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This paper discusses the application of robust control theory to evaluate 
robustness of the Ares-I control systems. Three techniques for estimating 
upper and lower bounds of uncertain parameters which yield stable closed- 
loop response are used here: (1) Monte Carlo analysis, (2) p analysis, and (3) 
characteristic frequency response analysis. All three methods are used to 
evaluate stability envelopes of the Ares-I control systems with uncertain 
aerodynamics and flex dynamics. The results show that characteristic 
frequency response analysis is the most effective of these methods for 
assessing robustness. 
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II. Introduction 


The Ares-I autopilot is a complex, time-varying control system. Mass properties and 
aerodynamic forces vary considerably over the launch timeline. In addition, propellant slosh 
varies during upper stage flight as propellant is consumed. The autopilot’s robustness in the 
presence of parameter uncertainty must be accurately evaluated to provide a precise stability 
envelope for the launch vehicle control systems. The stability envelope can be used as a 
yardstick for system engineers to determine if control system redesign is necessary for a given 
set of flight parameter variations such as aerodynamics and flex characteristics. 

The most common method for assessing control systems’ robustness is Monte Carlo 
analysis 11 , which is essentially a repeated guess and check screening process. Uncertainty in 
parameter values is modeled as a variation from nominal. Random samples of these parameters 
are taken to form unique system configurations or dispersions. The stability of each dispersion is 
determined through eigenanalysis, resulting in a map delineating stabilizing and destabilizing 
system configurations. 

An alternative way to handle uncertainty in system parameters is by developing an uncertain 
system model using robust control theory 2 ' 9 . Unlike the Monte Carlo approach, in which the 
system is re-formed for each set of parameter variations, the uncertain system model is formed 
only once. Wise considers several such techniques as applied to a simplified missile autopilot 
system in his paper on robustness tests 9 . In this paper, /j, analysis and characteristic frequency 
response analysis are used to estimate the perturbations that the system can tolerate before 
becoming unstable. An iterative // analysis algorithm is proposed for reducing conservatism 
when estimating bounds on uncertain parameters. Though certainly requiring fewer system 
formulations than Monte Carlo analysis, the algorithm’s iterative nature makes it costly in terms 
of computation and time. This provides motivation to develop an even more efficient method of 
analysis — characteristic frequency response analysis 6 . The characteristic frequency response 
method is the multi-input, multi-output (MIMO) version of Bode analysis. This algorithm 
represents a major improvement over other MIMO analysis methods such as Monte Carlo or // 
in that stability bounds can be determined using a single system formulation. 

In this paper, all three methods are used to evaluate stability envelopes of the Ares-I control 
systems with uncertain aerodynamics and flex dynamics. The iterative ju analysis technique 
offers a systematic way to detect stability bounds. Unlike Monte Carlo, which randomly selects 
parameter values, ju analysis uses variable step sizes in an optimized directional search. Since // 
requires that the system be evaluated at frequencies across the spectrum of interest, each iteration 
involves multiple system realizations. Depending on how quickly the n algorithm converges, it 
may require more or less time than comparable Monte Carlo analysis. In the Ares-I case study, /i 
analysis takes longer than Monte Carlo. The characteristic frequency response algorithm is 
significantly faster than ju or Monte Carlo. This is because no iteration or random sampling is 
required. The basic algorithm is the fastest of those considered here and is excellent for 
obtaining a “first blush” estimate on stability bounds. In some systems, this estimate will reflect 
the actual bounds. However, in problematic systems that exhibit coupling between uncertain 
parameters, the modified algorithm can be used to accurately assess bounds. This carries the 
added benefit of differentiating between bounds on each parameter. This unique aspect of the 
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characteristic frequency response algorithm makes it a valuable tool in system design and sets it 
apart from many other robustness tests. 

III. Monte Carlo Analysis 

The most common method for assessing robustness, Monte Carlo analysis is essentially a 
repeated guess and check analysis procedure. Uncertainty in parameter values is modeled as a 
variation from nominal. Random samples of these parameters are taken to form unique system 
configurations or dispersions. The stability of each dispersion is determined through 
eigenanalysis, resulting in a map delineating stabilizing and destabilizing system configurations. 

The first step in performing Monte Carlo analysis is to identify the uncertain parameters to be 
studied. All other system parameters are held fixed while uncertain parameters are allowed to 
vary over specified ranges. Though any distribution may be used, uniform distribution is 
generally preferred because its lack of bias helps ensure adequate representation of the extreme 
values in the range. In this paper, all distributions are presumed to be uniform. 

Each uncertain parameter forms an axis in the sample space. The boundaries of the sample 
space are defined by the range of variation for the respective parameter. Consider the spring 
mass damper system shown in Fig. la. Suppose there is some uncertainty in its stiffness and 
damping. These parameters form axes in the sample space. For demonstration purposes, 
stiffness and damping ranges are allowed to include negative values. The resulting sample space 
is depicted in Fig. lb. 
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a) b) 

Fig. 1 Spring Mass Damper (a) and Corresponding Sample Space (b) 


As the number of uncertain parameters increases, the process is repeated to form a hyper cubic 
sample space. The uniform distribution of each parameter leads to an unbiased joint distribution 
in the sample space. 

Next, parameter values are randomly selected from the sample space. A unique dispersion is 
formed for each sampling. Each dispersion represents a possible system configuration within the 
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specified parameter ranges. The stability of these dispersions is then assessed through 
eigenanalysis. Finally, unstable dispersions are mapped into the sample space by plotting points 
representing their unique parameter combinations. The resulting scatter plot is a map of unstable 
regions (and conversely, stable regions) of the sample space. Fig. 2 shows the results of Monte 
Carlo analysis on the spring mass damper system from Fig. 1. 
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Fig. 2 Monte Carlo Scatter Plot 


As expected, Fig. 2 shows that negative values of damping and stiffness cause instability. 

IV. Uncertain System Model 

An alternative way to handle uncertainty in system parameters is by developing an uncertain 
system model. Unlike the Monte Carlo approach, in which the system is re-formed for each set 
of parameter variations, the uncertain system model is formed only once. Then, one of several 
methods is used to estimate the perturbations that the system can tolerate before becoming 
unstable. Two of these methods, // analysis and characteristic frequency response, are outlined 
below. 

The basis of the uncertain system model lies in representing uncertainty as a perturbation 
from nominal. Consider a generic uncertain parameter a with nominal value ao. Variations from 
nominal are represented by 8 such that 


a = ao + 5 (1) 

where 8 can assume positive or negative values. This formulation applies to both scalar and 
matrix parameters. Returning to the previously described spring mass damper system from Fig. 
la, uncertainty in c and k can be represented as perturbations from nominal. The resulting 
system is depicted in Fig. 3a. This uncertain system can also be represented in block diagram 
form as shown in Fig. 3b. 
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Fig. 3 Uncertain Spring Mass Damper and Corresponding Block Diagram 


The process of representing uncertain parameters as perturbations from nominal is helpful in 
visualizing sources of uncertainty in the system. However, its use is limited because 
perturbations are embedded in the system structure. In order to be more useful for robustness 
testing, perturbations (and hence uncertainty) must be separated from the rest of the system. One 
way to accomplish this is by replacing uncertain block elements with equivalent structures 
known as linear fractional transformations (LFTs) 2 " 4 . 

Consider a generic uncertain block a representing Eq. (1) with input u and output y. The 
resulting expression for y is 

y = (a 0 +5)u. (2) 

The equivalent LFT representation of a is given by Fig. 4 
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Fig. 4 Linear Fractional Transformation (LFT) of a 


where A is a partitioned matrix and A a is a matrix or scalar (depending on the structure of the 
uncertain parameter). A new expression is found for y in terms of A and A a as follows 2 " 4 


y - (A 2 2 + MlAa (I - AnAa) ' A]2 ) U , 


( 3 ) 


where A a and the elements of A are selected such that Eq. (3) is equivalent to Eq. (2). 

This process is repeated for every uncertain parameter and the resulting perturbation blocks 
are collected in the block diagonal structure A such that A = diag(A\, A 2 , . . . , A n ). This separates 
the perturbations from the rest of the system, M, leading to the familiar block diagonal 
uncertainty formulation shown in Fig. 5. 
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Fig. 5 Block Diagonal Uncertain System Formulation 


The above process can be applied to systems of varying size and complexity. Though some 
effort is initially required to build the uncertain system model, once developed it has several 
applications, including // analysis and characteristic frequency response techniques for assessing 
robustness. These techniques are described in the following sections. 

V. // Analysis 

Having formulated the uncertain system model in such a way that uncertainty is isolated from 
the rest of the system, the next step is determining the range of allowable perturbations to 
maintain stability. Conversely, this process can be thought of as determining the smallest gain 
applied to A that drives the system unstable 2 . This is done using ju(M), the structured singular 
value (SSV) of M, which is defined as 

— — = inf {ct( A)s.t. det(7 - AM) = 0} (4) 

ju(M) AeA 

where A is the set of all possible A 5 . If no A satisfies det(7 - AM) = 0 , then i u(M) = 0. 

The first step in determining ju(M) is to evaluate M(ja>) at each frequency of interest in the 
spectrum over which the system operates. Next, 1 / ju(M ) , which corresponds to the smallest A 
that destabilizes the system, is found at each frequency. Several popular algorithms, most 
employing an optimization scheme involving a search for suitable perturbations, exist for this 
purpose 3 . 

Next, ju(M) is plotted as a function of frequency and the peak value (3 is determined. Using 
ft , a test developed by Doyle and modified by Safonov shows that for any stabilizing A , 

<t(A) < k M (5) 


where 


k 




1 


P 


( 6 ) 
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This test can provide an estimate of the maximum allowable perturbation A stab by applying the 
gain k M to the current perturbation A such that 

= K A • (?) 

Here, the form of the perturbation must be considered. Perturbations can take on any value, 
including complex values. In practice, many perturbations are strictly real. However, the SSY 
function may contain discontinuities when only real perturbations are considered 7 . This problem 
can be alleviated by considering complex perturbations. Though the SSV function obtained from 
complex perturbations is continuous, its results are generally conservative 7 . 

In order to reduce conservatism, uncertain systems are “re-centered” by shifting the nominal 
values of uncertain parameters up or down. Consider a scaling factor A applied to the nominal 
values of uncertain parameters such that 


a = A a 0 + 8 . (8) 

Using this formulation, the system can be re-centered by simply increasing or decreasing A . 
This re-centering can be visualized by plotting the system’s center in the .s-plane in terms of A 
such that A = 1 represents the original system, A > 1 shifts the system up, and A < 1 shifts the 
system down. Strictly real perturbations appear as lines on the real axis extending outward from 
system center. Complex perturbations are represented by circles around the system center, as 
shown in Fig. 6. 

The first step in re-centering is to set A 0 = 1 (corresponding to the original uncertain system) 
and perform ju analysis. The resulting k u , which determines how much A can vary before 

driving the system unstable, defines a stable region around the system center as shown in Fig. 6a. 
Next, a new system center A l is selected from within the stable region. If attempting to 
determine an upper bound, A 1 is placed on the real axis near the right edge of the stable region, 
whereas the opposite is true when evaluating the lower bound. This process can be formalized as 
follows 


At- — 


nti(i+7^ 


Hi) 


for k — 0 

for k = 1, 2, ..., n 


( 9 ) 


- 1 < 7 < 1 

where y is a step size that also determines the direction of re-centering (i.e. y > 0 for upper 
bound search and y < 0 for lower bound search). Then, the system is re-formed and re- 
evaluated. A new k is obtained and the process is repeated until k becomes negligible (see 
Fig. 6b). 
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a) b) 

Fig. 6 Re-centering of Uncertain System to Reduce Conservatism 


The system center X n at which k M is negligible defines the bound on uncertain parameters as 
follows 


bound = A n - 1 . (10) 

This ju analysis algorithm is very effective at reducing conservatism in estimates of bounds 
on uncertain parameters. However, it does not distinguish between parameters, i.e. it does not 
provide information about which parameter causes instability. Though certainly requiring fewer 
system formulations than Monte Carlo analysis, the algorithm’s iterative nature makes it costly 
in terms of computation and time. This provides motivation to develop an even more efficient 
method of analysis — characteristic frequency response. 

VI. Characteristic Frequency Response 

The characteristic frequency response method is the multi-input, multi-output (MIMO) 
version of Bode analysis. The basic algorithm represents a major improvement over other 
MIMO analysis methods such as Monte Carlo or // in that stability bounds can be determined 
using a single system formulation. In addition, a modified characteristic frequency response 
algorithm, outlined in the next section, expands on the basic algorithm to address certain 
problematic systems and to identify bounds for each uncertain parameter individually. 

A. Basic Characteristic Frequency Response Algorithm 

Like Bode analysis, characteristic frequency response is performed on the open loop transfer 
function G(ja > ) . The closed loop uncertain system shown in Fig. 7a is re-cast in the familiar 
feedback form in Fig. 7b, with G(jco) labeled accordingly. 
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Fig. 7 Closed Loop Uncertain System in Original (a) and Feedback (b) Forms 

The closed loop transfer function (CLTF) can also be written in terms of G(jof) as follows 

CLTF = G( '-' C0) . (11) 

1 -G(jco) 

Unlike Bode analysis, in which the overall system magnitude and phase are assessed at 
frequencies across the spectrum, characteristic frequency response calculates the magnitude and 
phase of each eigenvalue of the system at each frequency. The resulting eigenvalues are then 
ordered from most dominant to least dominant in magnitude. One way of simultaneously 
viewing both gain and phase of a particular response is in a Nyquist plot. This also provides the 
basis for stability analysis using the Nyquist criterion. 

In its simplest form, the Nyquist criterion shows that Nyquist plots of the CLTF that encircle 
the origin indicate instability 10 . By considering the CLTF in terms of the open loop transfer 
function (OLTF), stability can be assessed using G(jco ) . In the case of positive feedback, as 
shown in Eq. (11), G(jco) encirclements of +1 are equivalent to CLTF encirclements of the 
origin. Similarly, G(jco) encirclements of -1 are equivalent to CLTF encirclements of the origin 
in systems with negative feedback 10 . Referring to Fig. 7 and the description of the uncertain 
system model, positive feedback corresponds to positive perturbations while negative feedback 
represents negative perturbations. 

The information in a Nyquist plot and the Nyquist criterion itself can be translated into 
appropriate forms for use in a Nichols plot. This is advantageous because Nichols plots are more 
intuitive than Nyquist plots and help spread the data out to form a clear picture of system 
behavior. The first step is to determine the gain and phase of the encirclement points. The point 
+1 occurs at unity magnitude and 0±360& degrees of phase. Similarly, the point -1 occurs at 
unity magnitude and 180 ± 360/c degrees of phase. 

Next, consider a gain ka applied to A . This represents a scaling factor on the perturbation. 
The original OLTF G(ja > ) becomes k G G(jco ) as shown in Fig. 8. 
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Fig. 8 Open Loop Transfer Function with Perturbation Gain Ug 


Because the aforementioned encirclement points both have unity magnitude, the system must be 
stable if k G \G(jco)\ < 1 at the phase corresponding to the encirclement points. In order to find the 
minimum destabilizing gain ko.crit, this relationship is manipulated such that 


k 


G,crit 


1 

G(jco) 


( 12 ) 


Rather than plotting |G(y<y)| in the y-axis as in typical Nichols plots, kc,crit is plotted. Phase 
makes up the x-axis as usual. 

The points of interest are those where kc,cnt crosses 0 or 180 ( ± 360k) degrees of phase. 
Hence forth, these phase crossings will be referred to as upper bound crossings and lower bound 
crossings, respectively. The value of kc,cnt at an upper bound crossing indicates how much 
positive gain can be applied to the specified perturbation A before the system goes unstable. 
Conversely, kc, C rit at a lower bound crossing is the amount of negative gain that drives the system 
unstable. When the plot intersects more than one upper or lower bound crossing, the minimum 
values (smallest absolute values) are used. 

Using the previous example of a spring mass damper with uncertainty in stiffness and 
damping, eigenvalues of the open loop uncertain system are evaluated. In this case, there is only 
one significant response. Next, kc, C rit is plotted vs. phase as shown in Fig. 9. 
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Characteristic Frequency Response 



Phase (deg) 


Fig. 9 Characteristic Frequency Response of Uncertain Spring Mass Damper 


There is a lower bound crossing with kG.crit =1.0 and no upper bound crossing. Hence, the 
system is stable for stiffness and damping variations between -100% and +co. This matches the 
results of the Monte Carlo analysis and the intuitive result that the system remains stable as long 
as stiffness and damping are positive. 

B. Modified Characteristic Frequency Response Algorithm 

Though the basic characteristic frequency response algorithm is effective at determining 
upper and lower bounds on uncertain parameters, it does not distinguish between parameters. In 
other words, the bounds derived from the basic algorithm correspond to the minimum possible 
variation in any parameter that causes instability. Since kc, is applied to all perturbations equally, 
the basic algorithm effectively performs a linear search in the parameter space with all 
parameters increasing or decreasing together. This is shown in Fig. 10, which depicts the 
effective search pattern corresponding to Fig. 9. 

Parameter Space (% Variation) 



Fig. 10 Effective Search Pattern of Basic Algorithm 
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The basic algorithm forms the first step of the modified algorithm. Having established upper 
and lower bounds for the entire system, the next step is to find specific bounds for each 
parameter. This is accomplished by selecting a bound and setting all but one of the uncertain 
parameters to that value. The result is a new system with a single uncertain parameter. The 
basic algorithm is used to evaluate this system, resulting in bounds specific to the uncertain 
parameter being examined. This process is repeated for all uncertain parameters and at both the 
upper and lower bounds. 


VII. Case Study 

In order to assess the performance of the analysis methods outlined in the previous section, a 
case study is conducted using the Ares-I crew launch vehicle. Ares-I is a 2-stage rocket designed 
to carry the Orion crew exploration vehicle (CEY) to low-Earth orbit. Throughout its launch 
timeline, Ares-I experiences changes in both mass properties and environment. This case study 
employs the aforementioned analysis techniques to determine Ares-I’ s robustness to these 
changes. 

In this case study, vehicle pitch is studied. Because Ares-I is nearly symmetric about its 
longitudinal axis, the results obtained from the pitch axis are equally applicable to yaw. Launch 
vehicle dynamics, including propellant slosh and structural bending, are described by Frosch and 
Vallely 11 . A PID controller is used as a stabilizing autopilot system. The resulting nominal 
system is cast in state space form such that 

x = Ax + Bu 

(13) 

y = Cx + Du 

Next, uncertain system parameters are selected for study. Of the many parameters with some 
uncertainty, aerodynamic forces, C za and vehicle bending mode frequencies, com are among the 
most sensitive to perturbations. For modeling purposes, perturbations in com are considered. 
Perturbations in these parameters are considered in the analysis that follows. 

A. Monte Carlo Analysis 

Monte Carlo analysis is performed using the nominal system model with many sampled 
parameter values from the sample space. In this case, the axes of the sample space are C za and 
com , which range from nominal value -80% to +80% and -125% to +200%, respectively. The 
stability scatter plot, shown in Fig. 1 1, is formed by analyzing 10,000 dispersions. The resulting 
stability bounds are -100% to +170% in C za and -68% in com 2 (no upper bound exists). 
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Unstable Parameter Dispersions 



C Zffi (% Variation) 


Fig. 11 Monte Carlo Scatter Plot for Uncertainty in Aerodynamics and Bending 

Frequencies 


B. // Analysis 

The uncertain system model, which separates perturbations in C za and a>bi 2 from the rest of the 
system, is used when performing n analysis. Using complex perturbations, the system center is 
initially scaled by X = 1 and k M is obtained through // analysis. Then, X is decreased and the 

process is repeated until k fl < 0.01 . This occurs at X n =0.33, corresponding to a lower bound of 
-67%. The procedure is repeated, this time increasing X at each iteration to arrive at X n = 1.77. 
The upper bound is thus +177%. 

Comparison with the results obtained from Monte Carlo analysis reveals a problematic feature 
of the system. Unlike the spring mass damper used in previous examples, Ares-I exhibits cross- 
coupling between uncertain parameters. Careful examination of Fig. 1 1 shows that as a>b 2 
increases, the upper bound on C za increases. The effective search pattern for ju analysis is the 
same as for the basic characteristic frequency response algorithm, so both uncertain parameters 
increase and decrease together. The result is that the upper bound estimate produced by // 
analysis is overly optimistic. 

C. Characteristic Frequency Response 

As in /i analysis, characteristic frequency response techniques are applied to the uncertain 
system model. Magnitude and phase of the eigenvalues of the uncertain system are evaluated, 
then ordered by magnitude. Next, kG,cnt is plotted against phase for each of the dominant 
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eigenvalues (in this case, there are 2). Upper and lower bounds are determined by taking the 
minimum value of ko.cnt at upper bound crossings and lower bound crossings, as shown in Fig. 
12. This analysis results in upper and lower bound estimates of +177% and -68%, respectively. 


Characteristic Frequency Response 



Fig. 12 Characteristic Frequency Response of Uncertain Ares-I Model 


As with the upper bound obtained by // analysis, the upper bound predicted by the basic 
characteristic frequency response algorithm is overly optimistic. Again, this is due to the 
algorithm’s effective search pattern, which examines only the upper right and lower left comers 
of the parameter space. In order to explore all 4 comers of the parameter space, the modified 

algorithm is employed. This results in separate bounds of -100% to +171% in C za and -68% in 

2 

(Obi ■ 

D. Performance Summary 

Though the ultimate objective of all of the algorithms presented here is to identify bounds on 
stability, speed and computational demand are also useful in assessing the relative strengths and 
weaknesses of each algorithm. Here, computational demand is evaluated in terms of number of 
system formulations required. Table 1 summarizes the results of each algorithm’s analysis of 
Ares-I using a Windows Vista PC running dual core 1.86 GHz processors. 
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Method 

Number of 
Systems Formed 

Time 

(sec) 

Upper 

Bound 

Lower 

Bound 

Monte Carlo 

10,000 

265 

+170% 

-68% 

^Analysis 

57 

1,163 

+177% 

-67% 

Basic Char. 
Freq. Response 

1 

6 

+177% 

-68% 

Modified Char. 
Freq. Response 

5 

30 

+171% 

-68% 


Table 1 Summary of Analysis Results 


The basic characteristic frequency response algorithm is both the fastest at 6 seconds and the 
least computationally demanding, requiring a single system formulation. However, like fi 
analysis, it overestimates the upper bound by 7%. The modified characteristic frequency 
response algorithm requires more time and system formulations than the basic algorithm, but 
generates a more accurate upper bound. 

VIII. Conclusion 

In this paper, the characteristic frequency response algorithm and an iterative version of // 
analysis were described and tested. Their performance was assessed with Monte Carlo analysis 
as a benchmark test. The iterative / u analysis technique offers a systematic way to detect stability 
bounds. Unlike Monte Carlo, which randomly selects parameter values, // analysis uses variable 
step sizes in an optimized directional search. Since ft requires that the system be evaluated at 
frequencies across the spectrum of interest, each iteration involves multiple system realizations. 
Depending on how quickly the fi algorithm converges, it may require more or less time than 
comparable Monte Carlo analysis. In the above Ares-I case study, ft analysis takes longer than 
Monte Carlo. 

The characteristic frequency response algorithm is significantly faster than fi or Monte Carlo. 
This is because no iteration or random sampling is required. The basic algorithm is the fastest of 
those considered here and is excellent for obtaining a “first blush” estimate on stability bounds. 
In some systems, this estimate will reflect the actual bounds. However, in problematic systems 
that exhibit coupling between uncertain parameters, the modified algorithm can be used to 
accurately assess bounds. This carries the added benefit of differentiating between bounds on 
each parameter. This unique aspect of the characteristic frequency response algorithm makes it a 
valuable tool in system design and sets it apart from many other robustness tests. 
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